#include <fstream>
#include <sstream>
#include <string>
#include <iostream>
#include <vector>
#include <cmath>
#include "dbg.h"
#include "TMV.h"
#include "PowerSpectrum.h"
using std::abs;

//#define SN
//#define BO
#define CMB
#define LENS

#define LAMBDA
//#define LIMITW  // -1.05 < w < -0.95
#define SMOOTH
//#define ZPIVOT 0.4

//#define DOOVERALLMAX
//#define DOMARG
#define DOPLOT
//#define DOQUADMAX
//#define WWA
//#define NonFlat
//#define LSSOMDE
//#define ALTNOCMB

//#define FIXOM03 // only considers range 0.28 < Om < 0.32
//#define POSW  // Limit to w > -1 at all times.
//#define MAXOMDELSS 0.5

#ifdef MAXOMDELSS
#ifndef LSSOMDE
#define LSSOMDE
#endif
#endif

// Set NAME
#ifdef WWA

#ifdef NonFlat
#define NAME1 "wwa"
#else
#define NAME1 "wwa2"
#endif

#else
#ifdef CMB

#ifdef LAMBDA
#define NAME1 "lambda"
//#define NAME1 "lambda_pd"
#else
#define NAME1 "wmap3"
//#define NAME1 "cmb"
//#define NAME1 "cfrs"
#endif // LAMBDA

#elif defined LENS

#ifdef LAMBDA
#define NAME1 "lens_lambda"
#else
#define NAME1 "lens"
#endif // LAMBDA

#elif defined BO

#ifdef LAMBDA
#define NAME1 "bo_lambda"
#else 
#define NAME1 "bo"
#endif

#elif defined SN

#ifdef LAMBDA
#define NAME1 "sn_lambda"
#else
#define NAME1 "sn"
#endif // LAMBDA

#endif // CMB

#endif // WWA

#ifdef NonFlat
#define NAME "nf" + NAME1
#else
#define NAME NAME1
#endif

// Define CMBORLENS, ALLTHREE
#ifdef CMB
#define CMBORLENS
#ifdef LENS
#ifdef SN
#define ALLTHREE
#endif // SN
#endif // LENS
#else // !CMB
#ifdef LENS
#define CMBORLENS
#endif // LENS
#endif // CMB

#ifdef LAMBDA
#define WIsM1
#endif
#ifdef LIMITW
#define WIsM1
#endif

// For central, use 300, 8
// For outer, use 200, 12
#ifdef ALLTHREE
const int plotgridn=300;
const double smoothing_scale=5.; // in grid pixels
#else
const int plotgridn=200;
const double smoothing_scale=5.; // in grid pixels
#endif
const int smooth_iter=0;

#ifdef WWA
#ifdef NonFlat
const int allgridn=8;
#else
const int allgridn=12;
#endif
#else
#ifdef NonFlat
const int allgridn=12;
#else
const int allgridn=20;
#endif
#endif
const int singlegridn=100;

#ifdef WWA
//#ifdef LENS
//#define SUPPRESSHIGHWA
//#endif
#endif

#ifdef SN
extern "C" {
  void compute_sn1a_chi2_(const double& Om_m, const double& Om_de, const double& w, const double& wa, double& chisq);
}
#endif

using std::endl;
using std::vector;
using std::string;

std::ostream* dbgout=0;
bool XDEBUG=false;
//bool XDEBUG=true;

const string likefile = string("lnlike.") + NAME + ".dat";
const string debugfile = string("binlike.") + NAME + ".debug";

string om_s_smoutfile = "/home/mjarvis/plots/om_s.image";
string om_s_pgoutfile = "/home/mjarvis/plots/om_s.list";
string lam_w_smoutfile = "/home/mjarvis/plots/lam_w.image";
string lam_w_pgoutfile = "/home/mjarvis/plots/lam_w.list";
string lam_wp_smoutfile = "/home/mjarvis/plots/lam_wp.image";
string lam_wp_pgoutfile = "/home/mjarvis/plots/lam_wp.list";
string w_wa_smoutfile = "/home/mjarvis/plots/w_wa.image";
string w_wa_pgoutfile = "/home/mjarvis/plots/w_wa.list";
string wp_wa_smoutfile = "/home/mjarvis/plots/wp_wa.image";
string wp_wa_pgoutfile = "/home/mjarvis/plots/wp_wa.list";
string omdelss_w_smoutfile = "/home/mjarvis/plots/omdelss_w.image";
string omdelss_w_pgoutfile = "/home/mjarvis/plots/omdelss_w.list";

inline double SQR(double x) { return x*x; }

void QuadMax(size_t i1, size_t i2, const vector<double>&  f,
    double* xmax, double* fmax)
// Given f(i) = f[i] for i=i1..i2 and f(x) = dx^2 + ex + f
// Find the maximum value of f, with f(xmax) = fmax
{
  dbg<<"Start quadmax: \n";
  for(size_t i=i1;i<=i2;i++) dbg<<"f("<<i<<") = "<<f[i]<<endl; 
  size_t npts = i2-i1+1;
  Assert(npts>=3);

  tmv::Matrix<double> A(npts,3,0.);
  tmv::Vector<double> B(npts,0.);
  for(size_t i=i1,j=0;i<=i2;i++,j++) {
    A(j,0) = i*i;  A(j,1) = i;  A(j,2) = 1.;  
    B(j) = f[i];
  }
  tmv::Vector<double> C = B/A;
  Assert(C.size() == 3);

  *xmax = -C(1)/C(0)/2.;
  tmv::Vector<double> X(3,0.);
  X(0) = (*xmax)*(*xmax);  X(1) = *xmax;  X(2) = 1.;
  *fmax = X*C;
}

void QuadMax(size_t i1, size_t i2, size_t j1, size_t j2,
    const vector<vector<double> >& f,
    double* imax, double* jmax, double* fmax)
// Given f(i,j) = f[i][j] for i=i1..i2, j=j1..j2
// and f(i,j) = Ai^2 + Bj^2 + Cij + Di + Ej + F
// Find the maximum value of f, with f(imax,jmax) = fmax
{
  dbg<<"Start 2D quadmax: \n";
  size_t npts = (i2-i1)*(j2-j1);
  Assert(npts>=6);

  dbg<<"npts = "<<npts<<endl;
  tmv::Matrix<double> A(npts,6,0.);
  tmv::Vector<double> B(npts,0.);
  size_t n=0;
  for(size_t i=i1;i<i2;i++) for(size_t j=j1;j<j2;j++) {
    double Bn = f[i][j];
    if (Bn==0.) continue;
    A(n,0) = i*i;
    A(n,1) = j*j;
    A(n,2) = i*j;
    A(n,3) = i;
    A(n,4) = j;
    A(n,5) = 1.;
    B(n) = Bn;
    n++;
  }

  dbg<<"made A,B\n";
  tmv::Vector<double> C = B.subVector(0,n)/A.subMatrix(0,n,0,6);
  dbg<<"solved for C\n";
  Assert(C.size() == 6);
  dbg<<"C = "<<C<<endl;

  // dfdi = 2Ai + Cj + D == 0
  // dfdj = 2Bj + Ci + E == 0
  // Solve 2 equations simultaneously:
 
  tmv::Matrix<double> D(2,2,0.);
  tmv::Vector<double> E(2,0.);

  D(0,0) = 2*C[0];
  D(0,1) = D(1,0) =C[2];
  D(1,1) = 2*C[1];
  E(0) = -C[3];
  E(1) = -C[4];

  dbg<<"D = "<<D<<endl;
  dbg<<"E = "<<E<<endl;

  tmv::Vector<double> F = E/D;

  dbg<<"F = E/D = "<<F<<endl;

  *imax = F(0); *jmax = F(1);

  tmv::Vector<double> X(6,0.);
  X(0) = (*imax)*(*imax);
  X(1) = (*jmax)*(*jmax);
  X(2) = (*imax)*(*jmax);
  X(3) = (*imax);
  X(4) = (*jmax);
  X(5) = 1.;
  *fmax = X*C;
  dbg<<"fmax = "<<*fmax<<endl;
}

void QuadMax(size_t i1, size_t i2, size_t j1, size_t j2, size_t k1, size_t k2,
    size_t l1, size_t l2, size_t m1, size_t m2, size_t n1, size_t n2,
    const vector<vector<vector<vector<vector<vector<double> > > > > >&  f,
    double* imax, double* jmax, double* kmax, double* lmax, double* mmax, double* nmax,
    double* fmax)
// Given f(i,j,k,l,m,n) = f[i][j][k][l][m][n] for i=i1..i2, etc.
// and f(i,j,k,l,m,n) = Ai^2 + Bj^2 + Ck^2 + Dl^2 + Em^2 + Fn^2 +
//                      Gij + Hik + Iil + Jim + Kin +
//                      Ljk + Mjl + Njm + Ojn +
//                      Pkl + Qkm + Rkn +
//                      Slm + Tln +
//                      Umn +
//                      Vi + Wj + Xk + Yl + Zm + AAn + BB
// Find the maximum value of f, with f(umax,vmax,...,zmax) = fmax
{
  dbg<<"Start 6D quadmax: \n";
  size_t npts = (i2-i1)*(j2-j1)*(k2-k1)*(l2-l1)*(m2-m1)*(n2-n1);
  Assert(npts>=28);

  dbg<<"npts = "<<npts<<endl;
  tmv::Matrix<double> A(npts,28,0.);
  dbg<<"created A"<<endl;
  dbg<<"sizeof A = "<<28*npts*sizeof(double)/1024/1024<<" MB\n";
  tmv::Vector<double> B(npts,0.);
  dbg<<"created B"<<endl;
  dbg<<"sizeof B = "<<npts*sizeof(double)/1024/1024<<" MB\n";
  vector<size_t> i(7);
  i[6] = 1;
  size_t n=0;
  for(i[0]=i1;i[0]<i2;i[0]++) for(i[1]=j1;i[1]<j2;i[1]++) {
    for(i[2]=k1;i[2]<k2;i[2]++) for(i[3]=l1;i[3]<l2;i[3]++)
      for(i[4]=m1;i[4]<m2;i[4]++) for(i[5]=n1;i[5]<n2;i[5]++) {
	//dbg<<"i0,i1,i2,i3,i4,i5 = "<<i[0]<<","<<i[1]<<","<<i[2]<<","<<i[3]<<","<<i[4]<<","<<i[5]<<endl;
	//dbg<<"n = "<<n<<endl;
	double Bn = f[i[0]][i[1]][i[2]][i[3]][i[4]][i[5]];
	if (Bn==0.) continue;
	for(size_t j=0,m=0;j<=6;j++) for(size_t k=j;k<=6;k++,m++) {
	  Assert(n<A.colsize());
	  Assert(m<A.rowsize());
	  A(n,m) = i[j] * i[k];
	}
	//dbg<<"A.row(n) = "<<A.row(n)<<endl;
	B(n) = Bn;
	//dbg<<"B(n) = "<<B(n)<<endl;
	n++;
      }
  }

  dbg<<"made A,B\n";
  tmv::Vector<double> C = B.subVector(0,n)/A.subMatrix(0,n,0,28);
  dbg<<"solved for C\n";
  Assert(C.size() == 28);
  dbg<<"C = "<<C<<endl;

  // dfdi = 2Ai + Gj + Hk + Il + Jm + Kn + V == 0
  // likewise for 5 other derivatives
  // Solve 6 equations simultaneously:

  tmv::Matrix<double> D(6,6,0.);
  tmv::Vector<double> E(6,0.);

  for(size_t j=0,m=0;j<=6;j++) for(size_t k=j;k<=6;k++,m++) {
    if (j==6) continue; // Deriv of constant term is 0
    else if (k==6) // Deriv of linear term goes on right side
      E(j) = -C[m];
    else if (k==j) // Deriv of ^2 term
      D(j,k) = 2*C[m];
    else { // Deriv of ij type term
      D(j,k) = C[m];
      D(k,j) = C[m];
    }
  }

  dbg<<"D = "<<D<<endl;
  dbg<<"E = "<<E<<endl;

  tmv::Vector<double> F = E/D;

  dbg<<"F = E/D = "<<F<<endl;

  *imax = F(0); *jmax = F(1); *kmax = F(2); *lmax = F(3); *mmax = F(4); *nmax = F(5);

  tmv::Vector<double> X(28,0.);
  for(size_t j=0,m=0;j<=6;j++) for(size_t k=j;k<=6;k++,m++) {
    X(m) = (j<6 ? F(j) : 1.) * (k<6 ? F(k) : 1.);
  }
  *fmax = X*C;
  dbg<<"fmax = "<<*fmax<<endl;
}

double LinearFit(size_t i1, size_t i2, const vector<double>& f, 
    double fx)
{
  dbg<<"Start linearfit: \n";
  dbg<<"i1 = "<<i1<<", i2 = "<<i2<<endl;
  for(size_t i=i1;i<=i2;i++) dbg<<"f("<<i<<") = "<<f[i]<<endl; 
  dbg<<"Solve for f(x) = "<<fx<<endl;
  size_t npts = i2-i1+1;
  Assert(npts > 2);

  tmv::Matrix<double> A(npts,2);
  tmv::Vector<double> B(npts);
  for(size_t i=i1,j=0;i<=i2;i++,j++) {
    A(j,0) = i;  A(j,1) = 1.;
    B(j) = f[i];
  }
  tmv::Vector<double> C = B/A;
  Assert(C.size()==2);
  return (fx-C(1))/C(0);
}

double QuadFit(size_t i1, size_t i2, const vector<double>& f, 
    double fx, size_t iguess)
// Given f(i) = f[i] for i=i1..i2 and f(x) = dx^2 + ex + f
// Find x such that f(x) = fx
{
  dbg<<"Start quadfit: \n";
  dbg<<"i1 = "<<i1<<", i2 = "<<i2<<endl;
  for(size_t i=i1;i<=i2;i++) dbg<<"f("<<i<<") = "<<f[i]<<endl; 
  dbg<<"Solve for f(x) = "<<fx<<endl;
  size_t npts = i2-i1+1;
  Assert(npts > 3);

  tmv::Matrix<double> A(npts,3);
  tmv::Vector<double> B(npts);
  for(size_t i=i1,j=0;i<=i2;i++,j++) {
    A(j,0) = i*i;  A(j,1) = i;  A(j,2) = 1.;  
    B(j) = f[i];
  }
  tmv::Vector<double> C = B/A;
  Assert(C.size()==3);

  double d = C(1)*C(1) - 4.*C(0)*(C(2)-fx);
  if (d < 0) {
    dbg<<"Warning: unsolvable quadratic found in QuadFit\n";
    dbg<<"Doing linear approx.\n";
    if (std::abs(int(i1)-int(iguess)) < std::abs(int(i2)-int(iguess)))
      return LinearFit(i1,2*iguess-i1,f,fx);
    else
      return LinearFit(2*iguess-i2,i2,f,fx);
  }
  double x1 = ( -C(1) + sqrt(d) )/C(0)/2.;
  double x2 = ( -C(1) - sqrt(d) )/C(0)/2.;
  tmv::Vector<double> X1(3,0.);
  X1(0) = x1*x1; X1(1) = x1;  X1(2) = 1.;
  dbg<<"Solution1: f("<<x1<<") = "<<X1*C<<endl;
  tmv::Vector<double> X2(3,0.);
  X2(0) = x2*x2; X2(1) = x2;  X2(2) = 1.;
  dbg<<"Solution2: f("<<x2<<") = "<<X2*C<<endl;

  bool x1ok = (x1 >= i1 && x1 <= i2);
  bool x2ok = (x2 >= i1 && x2 <= i2);
  if (x1ok) {
    if (x2ok) {
      dbg<<"Warning: two solutions found within range:\n";
      dbg<<"x1 = "<<x1<<", x2 = "<<x2<<endl;
      if (std::abs(x1-iguess) < std::abs(x2-iguess)) {
	dbg<<"Using x1, since closer to guess: "<<iguess<<endl;
	return x1;
      } else {
	dbg<<"Using x2, since closer to guess: "<<iguess<<endl;
	return x2;
      }
    } else return x1;
  } else {
    if (x2ok) return x2;
    else {
      dbg<<"Warning: no solutions found within range:\n";
      dbg<<"x1 = "<<x1<<", x2 = "<<x2<<endl;
      dbg<<"Doing linear approx.\n";
      if (abs(int(i1)-int(iguess)) < abs(int(i2)-int(iguess)))
	return LinearFit(i1,2*iguess-i1,f,fx);
      else
	return LinearFit(2*iguess-i2,i2,f,fx);
    }
  }
}

double FitLower(const vector<double>& bins, double besti, double val,
    size_t gridn, double xmin, double xmax)
{
  size_t i=0;
  while (bins[i] < val) {i++; Assert(i<bins.size()); }
  if (i<4) return -999.;
  else if (int(bins.size())-int(floor(besti))<=4) return 999.;
  else {
    double xi = QuadFit(i-4,size_t(floor(besti))+4,bins,val,i);
    return (xi+0.5)/gridn * (xmax-xmin) + xmin;
  }
}

double FitUpper(const vector<double>& bins, double besti, double val,
    size_t gridn, double xmin, double xmax)
{
  size_t i=bins.size()-1;
  while (bins[i] < val) {Assert(i>0); i--; }
  if (int(floor(besti))<4) return -999.;
  else if (int(bins.size())-i<=4) return 999.;
  else {
    double xi = QuadFit(size_t(floor(besti))-4,i+4,bins,val,i);
    return (xi+0.5)/gridn * (xmax-xmin) + xmin;
  }
}

class BinArray1 {

  public:

    BinArray1(size_t _gridn, double _xmin, double _xmax) :
      gridn(_gridn), xmin(_xmin), xmax(_xmax), bins(_gridn,0.)
#ifdef ALTNOCMB
	, count(_gridn,0)
#endif
    {}

    void Add(double x, double val)
    {
      if (x > xmin && x < xmax) {
	size_t i = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	Assert(i < gridn);

	bins[i] += val;
#ifdef ALTNOCMB
	count[i]++;
#endif
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      for(size_t i=0;i<gridn;++i) 
	if (count[i]) bins[i] /= count[i];
    }
#endif

    double FindMax(double* xfit=0, double* xm1=0, double* xp1=0, double* xm2=0, double* xp2=0)
    {
      double besti=0;
      double maxlike=0.;
      for(size_t i=0;i<gridn;i++) {
	if (bins[i] > maxlike) {
	  maxlike = bins[i];
	  besti = i;
	  dbg<<"maxlike = "<<maxlike<<" at "<<i<<endl;
	}
      }

#ifdef DOQUADMAX
      int n = 10;
      // Refit based on quadratic peak:
      if (besti<n) {
	dbg<<"besti = "<<besti<<" in BinArray1::FindMax.\n";
	QuadMax(0,int(besti)+n,bins,&besti,&maxlike);
	dbg<<"extrapolated besti = "<<besti<<endl;
	myerror("xmin not low enough to include likelihood peak");
      } else if (besti>=gridn-n) {
	dbg<<"besti = "<<besti<<" in BinArray1::FindMax.\n";
        QuadMax(int(besti)-n,gridn-1,bins,&besti,&maxlike);
	dbg<<"extrapolated besti = "<<besti<<endl;
	myerror("xmax not high enough to include likelihood peak");
      }
      dbg<<"Best fit at i = "<<besti<<" = "<<maxlike<<endl;
      Assert(besti > 0.);
      double newbesti;
      QuadMax(int(besti)-n,int(besti)+n,bins,&newbesti,&maxlike);
      dbg<<"New bestfit at "<<newbesti<<" = "<<maxlike<<endl;
      if (newbesti < besti-n || newbesti > besti+n) {
	myerror("QuadMax extrapolated: change range or gridding for ths parameter");
      }
      besti = newbesti;
#endif

      if (xfit) {
	*xfit = (besti+0.5)/gridn * (xmax-xmin) + xmin;
        dbg<<"x = "<<*xfit<<endl;
      }

      if (xm1 || xp1) {
	double l1 = maxlike * exp(-1.0/2.);
	if (xm1) *xm1 = FitLower(bins,besti,l1,gridn,xmin,xmax);
	if (xp1) *xp1 = FitUpper(bins,besti,l1,gridn,xmin,xmax);
      }
      if (xm2 || xp2) {
	double l2 = maxlike * exp(-4.0/2.);
	if (xm2) *xm2 = FitLower(bins,besti,l2,gridn,xmin,xmax);
	if (xp2) *xp2 = FitUpper(bins,besti,l2,gridn,xmin,xmax);
      }
      return maxlike;
    }

  private:

    size_t gridn;
    double xmin, xmax;
    vector<double> bins;
#ifdef ALTNOCMB
    vector<int> count;
#endif
};

class BinArray2 {

  public:

    BinArray2(size_t _gridn, double _xmin, double _xmax, 
	double _ymin, double _ymax, bool _xlog=false, bool _ylog=false) :
      gridn(_gridn), xlog(_xlog), ylog(_ylog),
      xmin(xlog?log(_xmin):_xmin), xmax(xlog?log(_xmax):_xmax),
      ymin(ylog?log(_ymin):_ymin), ymax(ylog?log(_ymax):_ymax),
      bins(_gridn,vector<double>(_gridn,0.))
#ifdef ALTNOCMB
      , count(_gridn,vector<int>(_gridn,0)) 
#endif
      {}

    void Add(double _x, double _y, double val)
    {
      double x = xlog?log(_x):_x;
      double y = ylog?log(_y):_y;
      if (x > xmin && x < xmax && y > ymin && y < ymax ) {
	size_t i = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	size_t j = size_t(floor( (y-ymin)/(ymax-ymin) * gridn ));
	Assert(i < gridn);
	Assert(j < gridn);

	bins[i][j] += val;
#ifdef ALTNOCMB
	count[i][j]++;
#endif
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      xdbg<<"Start Average:\n";
      for(size_t i=0;i<gridn;++i) for(size_t j=0;j<gridn;++j) {
	xdbg<<"i,j = "<<i<<','<<j;
	double x = (i+0.5)/gridn * (xmax-xmin) + xmin;
	if (xlog) x = exp(x);
	double y = (j+0.5)/gridn * (ymax-ymin) + ymin;
	if (ylog) y = exp(y);
	xdbg<<"   ( "<<x<<" , "<<y<<" )\n";
	xdbg<<"count[i][j] = "<<count[i][j]<<endl;
	xdbg<<"bins[i][j] = "<<bins[i][j]<<" --> ";
	if (count[i][j])
	  bins[i][j] /= count[i][j];
	xdbg<<bins[i][j]<<endl;
      }
    }
#endif

    double FindMax(double* x=0, double* y=0) const
    {
      double besti=0,bestj=0;
      double maxlike=0.;
      for(size_t i=0;i<gridn;i++) for(size_t j=0; j<gridn;j++) {
      //for(int i=gridn-1;i>=0;i--) for(int j=gridn-1;j>=0;j--) {
	//}
	if (bins[i][j] > maxlike) {
	  maxlike = bins[i][j];
	  besti = i;
	  bestj = j;
	  dbg<<"maxlike = "<<maxlike<<" at "<<i<<','<<j<<endl;
	  dbg<<"x = "<<(i+0.5)/gridn * (xmax-xmin) + xmin<<endl;
	  dbg<<"y = "<<(j+0.5)/gridn * (ymax-ymin) + ymin<<endl;
	}
      }
      dbg<<"Best fit at "<<besti<<','<<bestj<<endl;

#ifdef DOQUADMAX
      // Refit based on quadratic peak:
      int n=30;
      if (besti<n) {
	dbg<<"besti = "<<besti<<" in BinArray2::FindMax.\n";
	myerror("umin not low enough to include likelihood peak");
      } else if (besti>=gridn-n) {
	dbg<<"besti = "<<besti<<" in BinArray2::FindMax.\n";
	myerror("umax not high enough to include likelihood peak");
      }
      if (bestj<n) {
	dbg<<"bestj = "<<bestj<<" in BinArray2::FindMax.\n";
	myerror("vmin not low enough to include likelihood peak");
      } else if (bestj>=gridn-n) {
	dbg<<"bestj = "<<bestj<<" in BinArray2::FindMax.\n";
	myerror("vmax not high enough to include likelihood peak");
      }
      QuadMax(int(besti)-n,int(besti)+n+1,int(bestj)-n,int(bestj)+n+1,
	  bins,&besti,&bestj,&maxlike);
#endif

      if (x) {
	*x = (besti+0.5)/gridn * (xmax-xmin) + xmin;
	if (xlog) *x = exp(*x);
      }
      if (y) {
	*y = (bestj+0.5)/gridn * (ymax-ymin) + ymin;
	if (ylog) *y = exp(*y);
      }
      return maxlike;
    }

    void MakePlot(const string& smoutfile, const string& pgoutfile)
    {
      string tag = "_";
#ifdef LENS
      tag += "wl";
#endif
#ifdef CMB
      tag += "cmb";
#endif
#ifdef SN
      tag += "sn";
#endif
#ifdef BO
      tag += "bo";
#endif

      double xmin_pg = xlog ? xmin/log(10.) : xmin;
      double xmax_pg = xlog ? xmax/log(10.) : xmax;
      double ymin_pg = ylog ? ymin/log(10.) : ymin;
      double ymax_pg = ylog ? ymax/log(10.) : ymax;

#ifdef SMOOTH
      // Smooth grid
      double maxval=0.;
      dbg<<"smoothing_scale = "<<smoothing_scale<<endl;
      const double smoothing_scale_sq = SQR(smoothing_scale);
      dbg<<"smoothing_scale_sq = "<<smoothing_scale_sq<<endl;
      int ismooth = int(3*smoothing_scale);
      dbg<<"ismooth = "<<ismooth<<endl;
      vector<vector<double> > newbins(gridn,vector<double>(gridn,0.));
      for(int i=0;i<int(gridn);i++) for(int j=0;j<int(gridn);j++) {
	xdbg<<"bins["<<i<<"]["<<j<<"] = "<<bins[i][j]<<endl;
	int i1 = std::max(i-ismooth,0);
	int i2 = std::min(i+ismooth+1,int(gridn));
	int j1 = std::max(j-ismooth,0);
	int j2 = std::min(j+ismooth+1,int(gridn));
	double newbin = 0.;
	double sumwt = 0.;
	for(int ii=i1;ii<i2;ii++) for(int jj=j1;jj<j2;jj++) {
	  //xdbg<<"ii,jj = "<<ii<<','<<jj<<endl;
	  double rr = SQR(ii-i) + SQR(jj-j);
	  rr /= smoothing_scale_sq;
	  //xdbg<<"rr = "<<rr<<endl;
	  double weight = exp(-rr/2.);
	  //xdbg<<"weight = "<<weight<<endl;
	  newbin += bins[ii][jj] * weight;
	  sumwt += weight;
	}
	newbin /= sumwt;
	if (newbin > maxval) maxval = newbin;
	newbins[i][j] = newbin;
	xdbg<<"newbins["<<i<<"]["<<j<<"] = "<<newbins[i][j]<<endl;
      }
      bins = newbins;
      dbg<<"Done gaussian convolution.\n";

      dbg<<"Now do inverse gradient smoothing\n";
      double maxwt = 2.e4/maxval; 
      // Essentially this rescales values so that maximum is 10^4
      // and differences of less than 1 are considered equal for the 
      // purposes of the weights.
      for(int iter=0;iter<smooth_iter;iter++) {
	vector<vector<double> > newbins(gridn,vector<double>(gridn,0.));
	for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++) {
	  size_t ii1= i==0 ? 0 : i-1;
	  size_t ii2= i==gridn-1 ? gridn-1 : i+1;
	  size_t jj1= j==0 ? 0 : j-1;
	  size_t jj2= j==gridn-1 ? gridn-1 : j+1;
	  double totwt=0.;
	  for(size_t ii=ii1;ii<=ii2;ii++) for(size_t jj=jj1;jj<=jj2;jj++) {
	    double weight = ( (i==ii && j==jj) || bins[i][j] == bins[ii][jj]) ?
	      maxwt : 1./abs(bins[i][j]-bins[ii][jj]);
	    if (!(weight <= maxwt)) weight = maxwt;
	    newbins[i][j] += bins[ii][jj] * weight;
	    totwt += weight;
	  }
	  if (totwt > 0.) newbins[i][j] /= totwt;
	  //xdbg<<"newbins["<<i<<"]["<<j<<"] = "<<newbins[i][j]<<"  totwt = "<<totwt/maxwt<<endl;
	}
	bins = newbins;
      }
      dbg<<"done\n";
#endif
	
      // Convert bins from likelihood to relative chisq
      double xmax,ymax;
      double maxlike = FindMax(&xmax,&ymax);
      double lnmaxlike = log(maxlike);
      dbg<<"Max at "<<xmax<<", "<<ymax<<" = "<<maxlike<<endl;
      for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++) {
	if (bins[i][j] > 0.1*maxlike) {
	  xdbg<<"bins["<<i<<"]["<<j<<"] = "<<bins[i][j]<<endl;
	}
	if (bins[i][j]) bins[i][j] = -2*log(bins[i][j]/maxlike);
	else bins[i][j] = 1.e10;
	//if (bins[i][j] < 4.605170186) {
	  //xdbg<<"chisq["<<i<<"]["<<j<<"] = "<<bins[i][j]<<endl;
	//}
      }

      // Output files for plotting programs

      std::ofstream smout((smoutfile+tag).c_str());
      smout.write((const char*)&gridn,4);  Assert(sizeof(gridn) == 4);
      smout.write((const char*)&gridn,4);  Assert(sizeof(gridn) == 4);
      //std::ofstream pgout((pgoutfile+tag).c_str());
      //pgout << xmin_pg<<' '<<xmax_pg<<' ';
      //pgout << ymin_pg<<' '<<ymax_pg<<endl;
      //pgout << gridn <<"  "<< gridn << endl;
      for(size_t n2=0;n2<gridn;n2++) for(size_t n1=0;n1<gridn;n1++) {
	xdbg<<"bins["<<n1<<"]["<<n2<<"] = "<<bins[n1][n2]<<endl;
	float val = bins[n1][n2];
	smout.write((const char*)&val,4); Assert(sizeof(val) == 4);
	//pgout << val <<endl;
      }
      smout.close();
      //pgout.close();
    }

  private:

    size_t gridn;
    bool xlog,ylog;
    double xmin, xmax;
    double ymin, ymax;
    vector<vector<double> > bins;
#ifdef ALTNOCMB
    vector<vector<int> > count;
#endif
};

class BinArray5 {

  public:

    BinArray5(size_t _gridn, double _umin, double _umax,
	double _vmin, double _vmax,
	double _wmin, double _wmax,
	double _xmin, double _xmax,
	double _ymin, double _ymax) :
      gridn(_gridn), umin(_umin), umax(_umax),
      vmin(_vmin), vmax(_vmax),
      wmin(_wmin), wmax(_wmax),
      xmin(_xmin), xmax(_xmax),
      ymin(_ymin), ymax(_ymax),
      bins(_gridn,vector<vector<vector<vector<double> > > >(_gridn,
	    vector<vector<vector<double> > >(_gridn,
	    vector<vector<double> >(_gridn,
	    vector<double>(_gridn,0.)))))
#ifdef ALTNOCMB
      , count(_gridn,vector<vector<vector<vector<int> > > >(_gridn,
	    vector<vector<vector<int> > >(_gridn,
	    vector<vector<int> >(_gridn,
	    vector<int>(_gridn,0))))) 
#endif
    {}

    void Add(double u, double v, double w, double x, double y, 
	double val)
    {
      if (u > umin && u < umax && v > vmin && v < vmax &&
	  w > wmin && w < wmax && x > xmin && x < xmax &&
	  y > ymin && y < ymax) {
	size_t i = size_t(floor( (u-umin)/(umax-umin) * gridn ));
	size_t j = size_t(floor( (v-vmin)/(vmax-vmin) * gridn ));
	size_t k = size_t(floor( (w-wmin)/(wmax-wmin) * gridn ));
	size_t l = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	size_t m = size_t(floor( (y-ymin)/(ymax-ymin) * gridn ));
	Assert(i < gridn);
	Assert(j < gridn);
	Assert(k < gridn);
	Assert(l < gridn);
	Assert(m < gridn);

	bins[i][j][k][l][m] += val;
#ifdef ALTNOCMB
	count[i][j][k][l][m]++;
#endif
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	  for(size_t m=0;m<gridn;m++) {
	    if (count[i][j][k][l][m])
	      bins[i][j][k][l][m] /= count[i][j][k][l][m];
	  }
    }
#endif

    double FindMax(double* u, double* v, double* w, double* x, double* y)
    {
      double besti=0,bestj=0,bestk=0,bestl=0,bestm=0,bestn=0;
      double maxlike=0.;
      for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	  for(size_t m=0;m<gridn;m++) {
	    if (bins[i][j][k][l][m] > maxlike) {
	      maxlike = bins[i][j][k][l][m];
	      besti = i;
	      bestj = j;
	      bestk = k;
	      bestl = l;
	      bestm = m;
	      dbg<<"maxlike = "<<maxlike<<" at "<<i<<','<<j<<','<<k<<','<<l<<','<<m<<endl;
	    }
	  }
      dbg<<"Best fit at "<<besti<<','<<bestj<<','<<bestk<<','<<bestl<<',';
      dbg<<bestm<<" = "<<maxlike<<endl;

      *u = (besti+0.5)/gridn * (umax-umin) + umin;
      *v = (bestj+0.5)/gridn * (vmax-vmin) + vmin;
      *w = (bestk+0.5)/gridn * (wmax-wmin) + wmin;
      *x = (bestl+0.5)/gridn * (xmax-xmin) + xmin;
      *y = (bestm+0.5)/gridn * (ymax-ymin) + ymin;

      return maxlike;
    }

  private:

    size_t gridn;
    double umin, umax;
    double vmin, vmax;
    double wmin, wmax;
    double xmin, xmax;
    double ymin, ymax;
    vector<vector<vector<vector<vector<double> > > > > bins;
#ifdef ALTNOCMB
    vector<vector<vector<vector<vector<int> > > > > count;
#endif
};

class BinArray6 {

  public:

    BinArray6(size_t _gridn, double _umin, double _umax,
	double _vmin, double _vmax,
	double _wmin, double _wmax,
	double _xmin, double _xmax,
	double _ymin, double _ymax,
	double _zmin, double _zmax) :
      gridn(_gridn), umin(_umin), umax(_umax),
      vmin(_vmin), vmax(_vmax),
      wmin(_wmin), wmax(_wmax),
      xmin(_xmin), xmax(_xmax),
      ymin(_ymin), ymax(_ymax),
      zmin(_zmin), zmax(_zmax),
      bins(_gridn,vector<vector<vector<vector<vector<double> > > > >(_gridn,
	    vector<vector<vector<vector<double> > > >(_gridn,
	    vector<vector<vector<double> > >(_gridn,
	    vector<vector<double> >(_gridn,
	    vector<double>(_gridn,0.))))))
#ifdef ALTNOCMB
      , count(_gridn,vector<vector<vector<vector<vector<int> > > > >(_gridn,
	    vector<vector<vector<vector<int> > > >(_gridn,
	    vector<vector<vector<int> > >(_gridn,
	    vector<vector<int> >(_gridn,
	    vector<int>(_gridn,0)))))) 
#endif
    {}

    void Add(double u, double v, double w, double x, double y, double z,
	double val)
    {
      if (u > umin && u < umax && v > vmin && v < vmax &&
	  w > wmin && w < wmax && x > xmin && x < xmax &&
	  y > ymin && y < ymax && z > zmin && z < zmax ) {
	size_t i = size_t(floor( (u-umin)/(umax-umin) * gridn ));
	size_t j = size_t(floor( (v-vmin)/(vmax-vmin) * gridn ));
	size_t k = size_t(floor( (w-wmin)/(wmax-wmin) * gridn ));
	size_t l = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	size_t m = size_t(floor( (y-ymin)/(ymax-ymin) * gridn ));
	size_t n = size_t(floor( (z-zmin)/(zmax-zmin) * gridn ));
	Assert(i < gridn);
	Assert(j < gridn);
	Assert(k < gridn);
	Assert(l < gridn);
	Assert(m < gridn);
	Assert(n < gridn);

	bins[i][j][k][l][m][n] += val;
#ifdef ALTNOCMB
	count[i][j][k][l][m][n]++;
#endif
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	  for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	    if (count[i][j][k][l][m][n])
	      bins[i][j][k][l][m][n] /= count[i][j][k][l][m][n];
	  }
    }
#endif

    double FindMax(double* u, double* v, double* w, double* x, double* y, double* z)
    {
      double besti=0,bestj=0,bestk=0,bestl=0,bestm=0,bestn=0;
      double maxlike=0.;
      for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	  for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	    if (bins[i][j][k][l][m][n] > maxlike) {
	      maxlike = bins[i][j][k][l][m][n];
	      besti = i;
	      bestj = j;
	      bestk = k;
	      bestl = l;
	      bestm = m;
	      bestn = n;
	      dbg<<"maxlike = "<<maxlike<<" at "<<i<<','<<j<<','<<k<<','<<l<<','<<m<<','<<n<<endl;
	    }
	  }
      dbg<<"Best fit at "<<besti<<','<<bestj<<','<<bestk<<','<<bestl<<',';
      dbg<<bestm<<','<<bestn<<" = "<<maxlike<<endl;

#ifdef XDOQUADMAX
      int nn=3;
      if (besti<nn) {
	dbg<<"besti = "<<besti<<" in BinArray6::FindMax.\n";
	myerror("umin not low enough to include likelihood peak");
      } else if (besti>=gridn-nn) {
	dbg<<"besti = "<<besti<<" in BinArray6::FindMax.\n";
	myerror("umax not high enough to include likelihood peak");
      }
      if (bestj<nn) {
	dbg<<"bestj = "<<bestj<<" in BinArray6::FindMax.\n";
	myerror("vmin not low enough to include likelihood peak");
      } else if (bestj>=gridn-nn) {
	dbg<<"bestj = "<<bestj<<" in BinArray6::FindMax.\n";
	myerror("vmax not high enough to include likelihood peak");
      }
      if (bestk<nn) {
	dbg<<"bestk = "<<bestk<<" in BinArray6::FindMax.\n";
	myerror("wmin not low enough to include likelihood peak");
      } else if (bestk>=gridn-nn) {
	dbg<<"bestk = "<<bestk<<" in BinArray6::FindMax.\n";
	myerror("wmax not high enough to include likelihood peak");
      }
      if (bestl<nn) {
	dbg<<"bestl = "<<bestl<<" in BinArray6::FindMax.\n";
	myerror("xmin not low enough to include likelihood peak");
      } else if (bestl>=gridn-nn) {
	dbg<<"bestl = "<<bestl<<" in BinArray6::FindMax.\n";
	myerror("xmax not high enough to include likelihood peak");
      }
      if (bestm<nn) {
	dbg<<"bestm = "<<bestm<<" in BinArray6::FindMax.\n";
	myerror("ymin not low enough to include likelihood peak");
      } else if (bestm>=gridn-nn) {
	dbg<<"bestm = "<<bestm<<" in BinArray6::FindMax.\n";
	myerror("ymax not high enough to include likelihood peak");
      }
      if (bestn<nn) {
	dbg<<"bestn = "<<bestn<<" in BinArray6::FindMax.\n";
	myerror("zmin not low enough to include likelihood peak");
      } else if (bestn>=gridn-nn) {
	dbg<<"bestn = "<<bestn<<" in BinArray6::FindMax.\n";
	myerror("zmax not high enough to include likelihood peak");
      }

      // Refit based on quadratic peak:
      QuadMax(0,gridn,0,gridn,0,gridn,0,gridn,0,gridn,0,gridn,
          bins,&besti,&bestj,&bestk,&bestl,&bestm,&bestn,&maxlike);
      dbg<<"New best fit at "<<besti<<','<<bestj<<','<<bestk<<','<<bestl<<',';
      dbg<<bestm<<','<<bestn<<" = "<<maxlike<<endl;
#endif

      *u = (besti+0.5)/gridn * (umax-umin) + umin;
      *v = (bestj+0.5)/gridn * (vmax-vmin) + vmin;
      *w = (bestk+0.5)/gridn * (wmax-wmin) + wmin;
      *x = (bestl+0.5)/gridn * (xmax-xmin) + xmin;
      *y = (bestm+0.5)/gridn * (ymax-ymin) + ymin;
      *z = (bestn+0.5)/gridn * (zmax-zmin) + zmin;

      return maxlike;
    }

  private:

    size_t gridn;
    double umin, umax;
    double vmin, vmax;
    double wmin, wmax;
    double xmin, xmax;
    double ymin, ymax;
    double zmin, zmax;
    vector<vector<vector<vector<vector<vector<double> > > > > > bins;
#ifdef ALTNOCMB
    vector<vector<vector<vector<vector<vector<int> > > > > > count;
#endif
};

class BinArray7 {

  public:

    BinArray7(size_t _gridn, 
	double _tmin, double _tmax,
	double _umin, double _umax,
	double _vmin, double _vmax,
	double _wmin, double _wmax,
	double _xmin, double _xmax,
	double _ymin, double _ymax,
	double _zmin, double _zmax) :
      gridn(_gridn),
      tmin(_tmin), tmax(_tmax),
      umin(_umin), umax(_umax),
      vmin(_vmin), vmax(_vmax),
      wmin(_wmin), wmax(_wmax),
      xmin(_xmin), xmax(_xmax),
      ymin(_ymin), ymax(_ymax),
      zmin(_zmin), zmax(_zmax),
      bins(_gridn,vector<vector<vector<vector<vector<vector<double> > > > > >(_gridn,
	    vector<vector<vector<vector<vector<double> > > > >(_gridn,
	    vector<vector<vector<vector<double> > > >(_gridn,
	    vector<vector<vector<double> > >(_gridn,
	    vector<vector<double> >(_gridn,
	    vector<double>(_gridn,0.)))))))
#ifdef ALTNOCMB
      , count(_gridn,vector<vector<vector<vector<vector<vector<int> > > > > >(_gridn,
	    vector<vector<vector<vector<vector<int> > > > >(_gridn,
	    vector<vector<vector<vector<int> > > >(_gridn,
	    vector<vector<vector<int> > >(_gridn,
	    vector<vector<int> >(_gridn,
	    vector<int>(_gridn,0))))))) 
#endif
    {}

    void Add(double t, double u, double v, double w, double x, double y, double z,
	double val)
    {
      if (t > tmin && t < tmax &&
	  u > umin && u < umax && v > vmin && v < vmax &&
	  w > wmin && w < wmax && x > xmin && x < xmax &&
	  y > ymin && y < ymax && z > zmin && z < zmax ) {
	size_t h = size_t(floor( (t-tmin)/(tmax-tmin) * gridn ));
	size_t i = size_t(floor( (u-umin)/(umax-umin) * gridn ));
	size_t j = size_t(floor( (v-vmin)/(vmax-vmin) * gridn ));
	size_t k = size_t(floor( (w-wmin)/(wmax-wmin) * gridn ));
	size_t l = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	size_t m = size_t(floor( (y-ymin)/(ymax-ymin) * gridn ));
	size_t n = size_t(floor( (z-zmin)/(zmax-zmin) * gridn ));
	Assert(h < gridn);
	Assert(i < gridn);
	Assert(j < gridn);
	Assert(k < gridn);
	Assert(l < gridn);
	Assert(m < gridn);
	Assert(n < gridn);

	bins[h][i][j][k][l][m][n] += val;
#ifdef ALTNOCMB
	count[h][i][j][k][l][m][n]++;
#endif
      } else {
	xdbg<<"Not in range: "<<t<<','<<u<<','<<v<<','<<w<<','<<x<<','<<y<<','<<z<<endl;
	if (!(t > tmin && t < tmax)) xdbg<<"t: "<<tmin<<" < "<<t<<" < "<<tmax<<" failed\n";
	if (!(u > umin && u < umax)) xdbg<<"u: "<<umin<<" < "<<u<<" < "<<umax<<" failed\n";
	if (!(v > vmin && v < vmax)) xdbg<<"v: "<<vmin<<" < "<<v<<" < "<<vmax<<" failed\n";
	if (!(w > wmin && w < wmax)) xdbg<<"w: "<<wmin<<" < "<<w<<" < "<<wmax<<" failed\n";
	if (!(x > xmin && x < xmax)) xdbg<<"x: "<<xmin<<" < "<<x<<" < "<<xmax<<" failed\n";
	if (!(y > ymin && y < ymax)) xdbg<<"y: "<<ymin<<" < "<<y<<" < "<<ymax<<" failed\n";
	if (!(z > zmin && z < zmax)) xdbg<<"z: "<<zmin<<" < "<<z<<" < "<<zmax<<" failed\n";
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      for(size_t h=0;h<gridn;h++)
	for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	  for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	    for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	      if (count[h][i][j][k][l][m][n])
		bins[h][i][j][k][l][m][n] /= count[h][i][j][k][l][m][n];
	    }
    }
#endif

    double FindMax(double* t, double* u, double* v, double* w, double* x, double* y, double* z)
    {
      double besth=0,besti=0,bestj=0,bestk=0,bestl=0,bestm=0,bestn=0;
      double maxlike=0.;
      for(size_t h=0;h<gridn;h++)
	for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	  for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	    for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	      if (bins[h][i][j][k][l][m][n] > maxlike) {
		maxlike = bins[h][i][j][k][l][m][n];
		besth = h;
		besti = i;
		bestj = j;
		bestk = k;
		bestl = l;
		bestm = m;
		bestn = n;
		dbg<<"maxlike = "<<maxlike<<" at "<<h<<','<<i<<','<<j<<','<<k<<','<<l<<','<<m<<','<<n<<endl;
	      }
	    }
      dbg<<"Best fit at ";
      dbg<<besth<<','<<besti<<','<<bestj<<','<<bestk<<','<<bestl<<',';
      dbg<<bestm<<','<<bestn<<" = "<<maxlike<<endl;

      *t = (besth+0.5)/gridn * (tmax-tmin) + tmin;
      *u = (besti+0.5)/gridn * (umax-umin) + umin;
      *v = (bestj+0.5)/gridn * (vmax-vmin) + vmin;
      *w = (bestk+0.5)/gridn * (wmax-wmin) + wmin;
      *x = (bestl+0.5)/gridn * (xmax-xmin) + xmin;
      *y = (bestm+0.5)/gridn * (ymax-ymin) + ymin;
      *z = (bestn+0.5)/gridn * (zmax-zmin) + zmin;

      return maxlike;
    }

  private:

    size_t gridn;
    double tmin, tmax;
    double umin, umax;
    double vmin, vmax;
    double wmin, wmax;
    double xmin, xmax;
    double ymin, ymax;
    double zmin, zmax;
    vector<vector<vector<vector<vector<vector<vector<double> > > > > > > bins;
#ifdef ALTNOCMB
    vector<vector<vector<vector<vector<vector<vector<int> > > > > > > count;
#endif
};

class BinArray8 {

  public:

    BinArray8(size_t _gridn, 
	double _smin, double _smax,
	double _tmin, double _tmax,
	double _umin, double _umax,
	double _vmin, double _vmax,
	double _wmin, double _wmax,
	double _xmin, double _xmax,
	double _ymin, double _ymax,
	double _zmin, double _zmax) :
      gridn(_gridn),
      smin(_smin), smax(_smax),
      tmin(_tmin), tmax(_tmax),
      umin(_umin), umax(_umax),
      vmin(_vmin), vmax(_vmax),
      wmin(_wmin), wmax(_wmax),
      xmin(_xmin), xmax(_xmax),
      ymin(_ymin), ymax(_ymax),
      zmin(_zmin), zmax(_zmax),
      bins(_gridn,
	  vector<vector<vector<vector<vector<vector<vector<double> > > > > > >(_gridn,
	    vector<vector<vector<vector<vector<vector<double> > > > > >(_gridn,
	      vector<vector<vector<vector<vector<double> > > > >(_gridn,
		vector<vector<vector<vector<double> > > >(_gridn,
		  vector<vector<vector<double> > >(_gridn,
		    vector<vector<double> >(_gridn,
		      vector<double>(_gridn,0.))))))))
#ifdef ALTNOCMB
      , count(_gridn,
	  vector<vector<vector<vector<vector<vector<vector<int> > > > > > >(_gridn,
	    vector<vector<vector<vector<vector<vector<int> > > > > >(_gridn,
	      vector<vector<vector<vector<vector<int> > > > >(_gridn,
		vector<vector<vector<vector<int> > > >(_gridn,
		  vector<vector<vector<int> > >(_gridn,
		    vector<vector<int> >(_gridn,
		      vector<int>(_gridn,0)))))))) 
#endif
    {}

    void Add(double s, double t, double u, double v, double w, 
	double x, double y, double z, double val)
    {
      if (s > smin && s < smax && t > tmin && t < tmax &&
	  u > umin && u < umax && v > vmin && v < vmax &&
	  w > wmin && w < wmax && x > xmin && x < xmax &&
	  y > ymin && y < ymax && z > zmin && z < zmax ) {
	size_t g = size_t(floor( (s-smin)/(smax-smin) * gridn ));
	size_t h = size_t(floor( (t-tmin)/(tmax-tmin) * gridn ));
	size_t i = size_t(floor( (u-umin)/(umax-umin) * gridn ));
	size_t j = size_t(floor( (v-vmin)/(vmax-vmin) * gridn ));
	size_t k = size_t(floor( (w-wmin)/(wmax-wmin) * gridn ));
	size_t l = size_t(floor( (x-xmin)/(xmax-xmin) * gridn ));
	size_t m = size_t(floor( (y-ymin)/(ymax-ymin) * gridn ));
	size_t n = size_t(floor( (z-zmin)/(zmax-zmin) * gridn ));
	Assert(g < gridn);
	Assert(h < gridn);
	Assert(i < gridn);
	Assert(j < gridn);
	Assert(k < gridn);
	Assert(l < gridn);
	Assert(m < gridn);
	Assert(n < gridn);

	bins[g][h][i][j][k][l][m][n] += val;
#ifdef ALTNOCMB
	count[g][h][i][j][k][l][m][n]++;
#endif
      } else {
	xdbg<<"Not in range: "<<s<<','<<t<<','<<u<<','<<v<<','<<w<<','<<x<<','<<y<<','<<z<<endl;
	if (!(s > smin && s < smax)) xdbg<<"s: "<<smin<<" < "<<s<<" < "<<smax<<" failed\n";
	if (!(t > tmin && t < tmax)) xdbg<<"t: "<<tmin<<" < "<<t<<" < "<<tmax<<" failed\n";
	if (!(u > umin && u < umax)) xdbg<<"u: "<<umin<<" < "<<u<<" < "<<umax<<" failed\n";
	if (!(v > vmin && v < vmax)) xdbg<<"v: "<<vmin<<" < "<<v<<" < "<<vmax<<" failed\n";
	if (!(w > wmin && w < wmax)) xdbg<<"w: "<<wmin<<" < "<<w<<" < "<<wmax<<" failed\n";
	if (!(x > xmin && x < xmax)) xdbg<<"x: "<<xmin<<" < "<<x<<" < "<<xmax<<" failed\n";
	if (!(y > ymin && y < ymax)) xdbg<<"y: "<<ymin<<" < "<<y<<" < "<<ymax<<" failed\n";
	if (!(z > zmin && z < zmax)) xdbg<<"z: "<<zmin<<" < "<<z<<" < "<<zmax<<" failed\n";
      }
    }

#ifdef ALTNOCMB
    void Average()
    {
      for(size_t g=0;g<gridn;g++) for(size_t h=0;h<gridn;h++)
	for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	  for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	    for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	      if (count[g][h][i][j][k][l][m][n])
		bins[g][h][i][j][k][l][m][n] /= count[g][h][i][j][k][l][m][n];
	    }
    }
#endif

    double FindMax(double* s, double* t, double* u, double* v, double* w, 
	double* x, double* y, double* z)
    {
      double bestg=0,besth=0,besti=0,bestj=0,bestk=0,bestl=0,bestm=0,bestn=0;
      double maxlike=0.;
      for(size_t g=0;g<gridn;g++) for(size_t h=0;h<gridn;h++)
	for(size_t i=0;i<gridn;i++) for(size_t j=0;j<gridn;j++)
	  for(size_t k=0;k<gridn;k++) for(size_t l=0;l<gridn;l++)
	    for(size_t m=0;m<gridn;m++) for(size_t n=0;n<gridn;n++) {
	      if (bins[g][h][i][j][k][l][m][n] > maxlike) {
		maxlike = bins[g][h][i][j][k][l][m][n];
		bestg = g;
		besth = h;
		besti = i;
		bestj = j;
		bestk = k;
		bestl = l;
		bestm = m;
		bestn = n;
		dbg<<"maxlike = "<<maxlike<<" at "<<g<<','<<h<<','<<i;
		dbg<<','<<j<<','<<k<<','<<l<<','<<m<<','<<n<<endl;
	      }
	    }
      dbg<<"Best fit at ";
      dbg<<bestg<<','<<besth<<','<<besti<<','<<bestj<<',';
      dbg<<bestk<<','<<bestl<<','<<bestm<<','<<bestn<<" = "<<maxlike<<endl;

      *s = (bestg+0.5)/gridn * (smax-smin) + smin;
      *t = (besth+0.5)/gridn * (tmax-tmin) + tmin;
      *u = (besti+0.5)/gridn * (umax-umin) + umin;
      *v = (bestj+0.5)/gridn * (vmax-vmin) + vmin;
      *w = (bestk+0.5)/gridn * (wmax-wmin) + wmin;
      *x = (bestl+0.5)/gridn * (xmax-xmin) + xmin;
      *y = (bestm+0.5)/gridn * (ymax-ymin) + ymin;
      *z = (bestn+0.5)/gridn * (zmax-zmin) + zmin;

      return maxlike;
    }

  private:

    size_t gridn;
    double smin, smax;
    double tmin, tmax;
    double umin, umax;
    double vmin, vmax;
    double wmin, wmax;
    double xmin, xmax;
    double ymin, ymax;
    double zmin, zmax;
    vector<vector<vector<vector<vector<vector<vector<vector<double> > > > > > > > bins;
#ifdef ALTNOCMB
    vector<vector<vector<vector<vector<vector<vector<vector<int> > > > > > > > count;
#endif
};


void ReadData(const string& likefile
#ifdef DOPLOT
    ,BinArray2& om_s_bins
#ifndef WIsM1
    ,BinArray2& lam_w_bins
#ifdef ZPIVOT
    ,BinArray2& lam_wp_bins
#endif
#endif
#ifdef WWA
    ,BinArray2& w_wa_bins
#ifdef ZPIVOT
    ,BinArray2& wp_wa_bins
#endif
#endif
#ifdef LSSOMDE
    ,BinArray2& omdelss_w_bins
#endif
#endif // DOPLOT
#ifdef DOOVERALLMAX
#ifdef WIsM1
    ,BinArray5& all_bins
#else
#ifdef WWA
#ifdef NonFlat
    ,BinArray8& all_bins
#else
    ,BinArray7& all_bins
#endif
#else
#ifdef NonFlat
    ,BinArray7& all_bins
#else
    ,BinArray6& all_bins
#endif
#endif // WWA
#endif // WIsM1
#endif // DOOVERALLMAX
#ifdef DOMARG
    ,BinArray1& omc_bins, BinArray1& omb_bins, BinArray1& h_bins,
    BinArray1& As_bins, BinArray1& n_bins,
#ifndef WIsM1
    BinArray1& w_bins,
#endif
#ifdef WWA
    BinArray1& wa_bins,
    BinArray1& wz_bins,
    //BinArray1& wzz_bins,
#endif
    BinArray1& Omm_bins, BinArray1& sig8_bins, BinArray1& sig8xx_bins
#endif // DOMARG
    )
{

  std::ifstream fin(likefile.c_str());
  if (!fin) myerror("Opening likefile ",likefile.c_str());

  string line;

  int linenum=0;
#ifdef WWA
  vector<size_t> wacount(15,0);
#endif
  while (getline(fin,line)) {
    if (linenum % 1000 == 0) {dbg<<linenum<<endl;}
    else {xdbg<<linenum<<"  "<<line<<endl;}
    linenum++;
    if (line[0] == '#') continue;

    double om_m,om_lam,om_b,h,As,n,w,wa=0.,sig8,datall,lnlike;

    std::istringstream linein(line);
#ifdef WWA
    linein >> om_m>>om_lam>>om_b>>h>>As>>n>>w>>wa>>sig8>>datall>>lnlike;
#else
    linein >> om_m>>om_lam>>om_b>>h>>As>>n>>w>>sig8>>datall>>lnlike;
#endif
    if (!linein) myerror("Reading lnlike.dat line",line.c_str());

#ifdef LIMITW
    if (w > -0.9 || w < -1.1) continue;
#endif
#ifdef FIXOM03
    if (om_m < 0.28 || om_m > 0.32) continue;
#endif
#ifdef POSW
    if (w < -1 || w+wa < -1) continue;
#endif

#ifdef LAMBDA
    w = -1.;
#endif

#ifdef LSSOMDE
    double alss = 1./(1.+1089.);
    double om_k = 1.-om_m-om_lam;
    double lna = log(alss);
    double ommlss1 = om_m*exp(-3.*lna);
    double omdelss1 = om_lam*exp(-3.*lna*(1+w+wa)-3.*wa*(1.-alss));
    double omklssl = om_k*exp(-2.*lna);
    double omtotlss1 = ommlss1+omdelss1+omklssl;
    double omdelss = omdelss1/omtotlss1;
    xdbg<<w<<"  "<<wa<<"  "<<om_lam<<"  "<<ommlss1<<","<<omdelss1<<"  "<<omdelss<<endl;
#ifdef MAXOMDELSS
    if (omdelss > MAXOMDELSS) continue;
#endif
#endif // LSSOMDE

    xdbg<<"Read lnlike = "<<lnlike<<", datall = "<<datall<<endl;

#ifdef ALTNOCMB
    lnlike = datall;
#endif

#ifdef CMB
#ifndef LENS
    lnlike -= datall;
    xdbg<<"after lens subtract: "<<lnlike<<endl;
#else
    // like *= exp(-(ln(om_m)-ln(0.5))^2/ln(2)^2)
    //lnlike -= 4.16*std::abs(log(om_m)+0.69);
#endif //!LENS
#endif //CMB

    xdbg<<"before sn\n";
#ifdef SN
#ifdef CMBORLENS
    double snchisq;
    xdbg<<"om_m = "<<om_m<<", om_de = "<<om_lam<<", w = "<<w<<endl;
    compute_sn1a_chi2_(om_m,om_lam,w,wa,snchisq);
    xdbg<<"SN chisq = "<<snchisq<<": lnlike = "<<lnlike<<" -> ";
    lnlike -= snchisq/2.;
    xdbg<<lnlike<<endl;
#endif // CMBORLENS
#endif // SN
    xdbg<<"after sn\n";

    xdbg<<"before bo\n";
#ifdef BO
    // From Eisenstein et al:
    // Baryon oscillations primarily constrain the ratio
    // R(z) = Dv(z)/Da(1089)
    // where Da(z) is the angular diameter distance
    // and Dv(z) = (Da(z)^2 cz/H(z))^1/3
    // for the z at which most of the galaxies lie.
    // The SDSS measurements are R(0.35) = 0.0979 +- 0.0036 (1-sigma)

    Universe u(om_m,om_lam,om_b,h,As,n,w,wa);
    xdbg<<"Made universe\n";
    double da1089 = u.fk(u.w_z(1089.));
    double da035 = u.fk(u.w_z(0.35));
    double h035 = u.hsq(1/1.35);
    Assert(h035 > 0.);
    h035 = sqrt(h035);
    double dv035 = da035*da035*0.35/h035;
    dv035 = exp(log(dv035)/3.);
    double r = dv035 / da1089;
    xdbg<<"for baryons: dv(0.35) = "<<dv035<<", da(1089) = "<<da1089<<endl;
    xdbg<<"ratio = "<<r<<": compare to 0.0979 += 0.0036\n";
    double bochisq = SQR(r-0.0979)/SQR(0.0036);
    xdbg<<"BO chisq = "<<bochisq<<endl;
    lnlike -= bochisq/2.;
    xdbg<<"new lnlike = "<<lnlike<<endl;
#endif
    xdbg<<"after bo\n";

    double like = exp(lnlike);
    xdbg<<"like = "<<like<<endl;

#ifdef DOPLOT
    xdbg<<"om_m = "<<om_m<<", sig8 = "<<sig8<<endl;
    om_s_bins.Add(om_m,sig8,like);
#ifndef WIsM1
    lam_w_bins.Add(om_lam,w,like);
#endif
#ifdef WWA
    //dbg<<"Wa: like("<<w<<","<<wa<<") = "<<like<<endl;
    w_wa_bins.Add(w,wa,like);
    wacount[int(floor(wa+5))]++;
#endif
#ifdef LSSOMDE
    omdelss_w_bins.Add(omdelss,w,like);
#endif

#ifdef ZPIVOT
    double wp = w + wa*(1.-1./(1.+ZPIVOT));
    lam_wp_bins.Add(om_lam,wp,like);
#ifdef WWA
    wp_wa_bins.Add(wp,wa,like);
#endif
#endif // ZPIVOT

#endif // DOPLOT


#ifdef DOOVERALLMAX
    all_bins.Add(om_m,
#ifdef NonFlat
	om_lam,
#endif
	om_b,h,As,n,
#ifndef WIsM1
	w,
#endif
#ifdef WWA
	wa,
#endif
	like);
#endif

#ifdef DOMARG
#ifdef SUPPRESSHIGHWA
    if (wa > 2.7) like *= 0.5;
#endif
    double hsq = h*h;
    omc_bins.Add((om_m-om_b)*hsq,like);
    omb_bins.Add(om_b*hsq,like);
    h_bins.Add(h,like);
    As_bins.Add(As,like);
    n_bins.Add(n,like);
#ifndef WIsM1
    w_bins.Add(w,like);
#endif
#ifdef WWA
    wa_bins.Add(wa,like);
#ifdef ZPIVOT
    wz_bins.Add(w+wa*(1.-1./(1.+ZPIVOT)),like);
    //wzz_bins.Add(w+wa*(1.-1./(1.+1.6)),like);
#endif
#endif
    Omm_bins.Add(om_m,like);
    sig8_bins.Add(sig8,like);
    sig8xx_bins.Add(sig8*exp(0.61*log(om_m/0.3)),like);
#endif
  }
  dbg<<"Done read: "<<linenum<<" lines\n";

#ifdef ALTNOCMB
#ifdef DOPLOT
  om_s_bins.Average();
#ifndef WIsM1
  lam_w_bins.Average();
#endif
#ifdef WWA
  w_wa_bins.Average();
#endif
#ifdef LSSOMDE
  omdelss_w_bins.Average();
#endif

#ifdef ZPIVOT
  lam_wp_bins.Average();
#ifdef WWA
  wp_wa_bins.Average();
#endif
#endif // ZPIVOT

#endif // DOPLOT

#ifdef DOOVERALLMAX
  all_bins.Average();
#endif

#ifdef DOMARG
  omc_bins.Average();
  omb_bins.Average();
  h_bins.Average();
  As_bins.Average();
  n_bins.Average();
#ifndef WIsM1
  w_bins.Average();
#endif
#ifdef WWA
  wa_bins.Average();
#ifdef ZPIVOT
  wz_bins.Average();
  //wzz_bins.Average();
#endif
#endif
  Omm_bins.Average();
  sig8_bins.Average();
  sig8xx_bins.Average();
#endif
#endif

#ifdef DOPLOT
#ifdef WWA
  for(size_t i=0;i<wacount.size();i++) dbg<<"wacount["<<i<<"] = "<<wacount[i]<<endl;
#endif
#endif

  if (linenum == 0) myerror("No lines read from ",likefile.c_str());
}

int main(int argc, char **argv)
{
  dbgout = new std::ofstream(debugfile.c_str());
  //dbgout = &std::cout;
  dbgout->setf(std::ios_base::unitbuf);

#ifdef DOPLOT
  BinArray2 om_s_bins(plotgridn, 0.1,0.6,0.45,1.5, true,true);
#ifndef WIsM1
#ifdef MAXOMDELSS
  BinArray2 lam_w_bins(plotgridn, 0.35,1.0,-2.0,0.0);
#else
  BinArray2 lam_w_bins(plotgridn, 0.35,1.0,-2.0,0.0);
#endif
#ifdef ZPIVOT
#ifdef MAXOMDELSS
  BinArray2 lam_wp_bins(plotgridn, 0.0,1.0,-1.0,0.0);
#else
  BinArray2 lam_wp_bins(plotgridn, 0.35,1.0,-1.0,0.0);
#endif
#endif
#endif // WIsM1
#ifdef WWA
#ifdef POSW
  BinArray2 w_wa_bins(plotgridn, -1.0,0.0,-1.0,1.0);
#else
#ifdef MAXOMDELSS
  BinArray2 w_wa_bins(plotgridn, -3.0,0.0,-2.0,4.0);
#else
  BinArray2 w_wa_bins(plotgridn, -3.0,0.0,-2.0,4.0);
#endif
#endif
#ifdef ZPIVOT
#ifdef MAXOMDELSS
  BinArray2 wp_wa_bins(plotgridn, -1.0,0.0,-2.0,6.0);
#else
  BinArray2 wp_wa_bins(plotgridn, -1.0,0.0,-2.0,4.0);
#endif
#endif
#endif // WWA

#ifdef LSSOMDE
#ifdef POSW
  BinArray2 omdelss_w_bins(plotgridn, 0.,0.15,-1.,0.);
#else
  BinArray2 omdelss_w_bins(plotgridn, 0.,1.0,-6.,0.);
#endif
#endif // LSSOMDE

#endif // DOPLOT

#ifdef DOOVERALLMAX
  //                             Om_m       Om_b          H       A_s       N          W
#ifdef WIsM1
  BinArray5 all_bins(allgridn, 0.2,0.3, 0.04,0.05, 0.65,0.75, 0.6,0.7, 0.9,1.0);
#else
#ifdef WWA
#ifdef NonFlat
  BinArray8 all_bins(allgridn, 0.27,0.35, 0.62,0.70, 0.046,0.054, 0.64,0.72, 0.75,0.83, 0.95,0.99, -1.3,-1.14, 1.14, 1.3 );
#else
  BinArray7 all_bins(allgridn, 0.35,0.45, 0.055,0.065, 0.60,0.70, 0.7,0.8, 0.95,1.00, -1.8,-1.3, 1.5, 2.0 );
#endif
#else
#ifdef NonFlat
  BinArray7 all_bins(allgridn, 0.2,0.3, 0.04,0.05, 0.65,0.75, 0.6,0.8, 0.9,1.0, -1.1,-0.8);
#else
  BinArray6 all_bins(allgridn, 0.2,0.3, 0.04,0.05, 0.65,0.75, 0.6,0.8, 0.9,1.0, -1.1,-0.8);
#endif
#endif
#endif
#endif

#ifdef DOMARG
  BinArray1 omc_bins(singlegridn, 0.0,0.3);
  BinArray1 omb_bins(singlegridn, 0.0,0.05);
  BinArray1 As_bins(singlegridn, 0.0,2.0);
  BinArray1 n_bins(singlegridn, 0.75,1.25);
#ifdef WWA
  BinArray1 h_bins(singlegridn, 0.3,1.1);
  BinArray1 Omm_bins(singlegridn, 0.1,0.9);
#ifdef LENS
  BinArray1 w_bins(singlegridn, -4.,0.);
  BinArray1 wa_bins(singlegridn, -2.,6.);
  BinArray1 wz_bins(singlegridn, -2.,0.);
  //BinArray1 wzz_bins(singlegridn, -2.,0.);
  BinArray1 sig_bins(singlegridn, 0.4,1.4);
  BinArray1 sigxx_bins(singlegridn, 0.4,1.4);
#else
  BinArray1 w_bins(singlegridn, -6.,0.);
  BinArray1 wa_bins(singlegridn, -2.,9.);
  BinArray1 wz_bins(singlegridn, -4.,0.);
  //BinArray1 wzz_bins(singlegridn, -1.5,1.);
  BinArray1 sig_bins(singlegridn, 0.0,1.4);
  BinArray1 sigxx_bins(singlegridn, 0.0,1.4);
#endif // LENS
#else // !WWA:
  BinArray1 h_bins(singlegridn, 0.5,0.9);
#ifndef WIsM1
  BinArray1 w_bins(singlegridn, -1.4,-0.6);
#endif
  BinArray1 Omm_bins(singlegridn, 0.1,0.5);
  BinArray1 sig_bins(singlegridn, 0.4,1.4);
  BinArray1 sigxx_bins(singlegridn, 0.4,1.4);
#endif // WWA
#endif // DOMARG

  dbg<<"before read "<<likefile<<endl;
  ReadData(likefile
#ifdef DOPLOT
      ,om_s_bins
#ifndef WIsM1
      ,lam_w_bins
#ifdef ZPIVOT
      ,lam_wp_bins
#endif
#endif
#ifdef WWA
      ,w_wa_bins
#ifdef ZPIVOT
      ,wp_wa_bins
#endif
#endif
#ifdef LSSOMDE
      ,omdelss_w_bins
#endif
#endif
#ifdef DOOVERALLMAX
      ,all_bins
#endif
#ifdef DOMARG
      ,omc_bins,omb_bins,h_bins,As_bins,n_bins,
#ifndef WIsM1
      w_bins,
#endif
#ifdef WWA
      wa_bins,
      wz_bins,
      //wzz_bins,
#endif
      Omm_bins,sig_bins,sigxx_bins
#endif
      );
  dbg<<"after read\n";

#ifdef DOOVERALLMAX
  {
    double bestomm,bestomb,besth,bestAs,bestn;
#ifdef NonFlat
    double bestomlam;
#endif
#ifndef WIsM1
    double bestw;
#endif
#ifdef WWA
    double bestwa;
#endif
    all_bins.FindMax(&bestomm,
#ifdef NonFlat
	&bestomlam,
#endif
	&bestomb,&besth,&bestAs,&bestn
#ifndef WIsM1
	,&bestw
#endif
#ifdef WWA
	,&bestwa
#endif
	);
    dbg<<"Overall best fit at:\n";
    dbg<<"Omega_m = "<<bestomm<<endl;
#ifdef NonFlat
    dbg<<"Omega_DE = "<<bestomlam<<endl;
#endif
    dbg<<"Omega_b = "<<bestomb<<endl;
    dbg<<"h = "<<besth<<endl;
    dbg<<"A_s = "<<bestAs<<endl;
    dbg<<"n = "<<bestn<<endl;
#ifndef WIsM1
    dbg<<"w = "<<bestw<<endl;
#endif
#ifdef WWA
    dbg<<"wa = "<<bestwa<<endl;
#endif
  }
#endif // DOOVERALLMAX

#ifdef DOPLOT
  dbg<<"Start make plot for omega-sig\n";
  om_s_bins.MakePlot(om_s_smoutfile,om_s_pgoutfile);
#ifndef WIsM1
  dbg<<"Start make plot for lambda-w\n";
  lam_w_bins.MakePlot(lam_w_smoutfile,lam_w_pgoutfile);
#ifdef ZPIVOT
  dbg<<"Start make plot for lambda-wp\n";
  lam_wp_bins.MakePlot(lam_wp_smoutfile,lam_wp_pgoutfile);
#endif
#endif
#ifdef WWA
  dbg<<"Start make plot for w-wa\n";
  w_wa_bins.MakePlot(w_wa_smoutfile,w_wa_pgoutfile);
#ifdef ZPIVOT
  dbg<<"Start make plot for wp-wa\n";
  wp_wa_bins.MakePlot(wp_wa_smoutfile,wp_wa_pgoutfile);
#endif
#endif
#ifdef LSSOMDE
  dbg<<"Start make plot for omdelss-w\n";
  omdelss_w_bins.MakePlot(omdelss_w_smoutfile,omdelss_w_pgoutfile);
#endif
#endif

#ifdef DOMARG
  double omc_m2,omc_m1,omc_fit,omc_p1,omc_p2;
  double omb_m2,omb_m1,omb_fit,omb_p1,omb_p2;
  double h_m2,h_m1,h_fit,h_p1,h_p2;
  double As_m2,As_m1,As_fit,As_p1,As_p2;
  double n_m2,n_m1,n_fit,n_p1,n_p2;
  double w_m2,w_m1,w_fit,w_p1,w_p2;
  double wa_m2,wa_m1,wa_fit,wa_p1,wa_p2;
  double wz_m2,wz_m1,wz_fit,wz_p1,wz_p2;
  //double wzz_m2,wzz_m1,wzz_fit,wzz_p1,wzz_p2;
  double Omm_m2,Omm_m1,Omm_fit,Omm_p1,Omm_p2;
  double sig_m2,sig_m1,sig_fit,sig_p1,sig_p2;
  double sigxx_m2,sigxx_m1,sigxx_fit,sigxx_p1,sigxx_p2;
  dbg<<"Find errors for omc\n";
  omc_bins.FindMax(&omc_fit,&omc_m1,&omc_p1,&omc_m2,&omc_p2);
  dbg<<"Find errors for omb\n";
  omb_bins.FindMax(&omb_fit,&omb_m1,&omb_p1,&omb_m2,&omb_p2);
  dbg<<"Find errors for h\n";
  h_bins.FindMax(&h_fit,&h_m1,&h_p1,&h_m2,&h_p2);
  dbg<<"Find errors for As\n";
  As_bins.FindMax(&As_fit,&As_m1,&As_p1,&As_m2,&As_p2);
  dbg<<"Find errors for n\n";
  n_bins.FindMax(&n_fit,&n_m1,&n_p1,&n_m2,&n_p2);
#ifndef WIsM1
  dbg<<"Find errors for w\n";
  w_bins.FindMax(&w_fit,&w_m1,&w_p1,&w_m2,&w_p2);
#endif
#ifdef WWA
  dbg<<"Find errors for wa\n";
  wa_bins.FindMax(&wa_fit,&wa_m1,&wa_p1,&wa_m2,&wa_p2);
#ifdef ZPIVOT
  dbg<<"Find errors for w(z="<<ZPIVOT<<")\n";
  wz_bins.FindMax(&wz_fit,&wz_m1,&wz_p1,&wz_m2,&wz_p2);
  //dbg<<"Find errors for w(z=1.6)\n";
  //wzz_bins.FindMax(&wzz_fit,&wzz_m1,&wzz_p1,&wzz_m2,&wzz_p2);
#endif
#endif
  dbg<<"Find errors for om\n";
  Omm_bins.FindMax(&Omm_fit,&Omm_m1,&Omm_p1,&Omm_m2,&Omm_p2);
  dbg<<"Find errors for sig\n";
  sig_bins.FindMax(&sig_fit,&sig_m1,&sig_p1,&sig_m2,&sig_p2);
  dbg<<"Find errors for sig(om/0.3)^1/2\n";
  sigxx_bins.FindMax(&sigxx_fit,&sigxx_m1,&sigxx_p1,&sigxx_m2,&sigxx_p2);

#ifdef WIsM1
  const double omb_m2x = 0.0191319;
  const double omb_p2x = 0.0254743;
  const double omc_m2x = 0.0964465;
  const double omc_p2x = 0.121955;
  const double n_m2x = 0.900396;
  const double n_p2x = 1.00446;
  const double As_m2x = 0.565400;
  const double As_p2x = 0.882033;
  const double h_m2x = 0.664930;
  const double h_p2x = 0.770375;
  const double Omm_m2x = 0.209035;
  const double Omm_p2x = 0.303747;
  const double sig_m2x = 0.725189;
  const double sig_p2x = 0.919465;
  const double sigxx_m2x = 0.;
  const double sigxx_p2x = 0.;
#else
#ifdef WWA
  /*
  const double omb_m2x = 0.0205416;
  const double omb_p2x = 0.0265663;
  const double omc_m2x = 0.0906327;
  const double omc_p2x = 0.130972;
  const double n_m2x = 0.922315;
  const double n_p2x = 1.059460;
  const double w_m2x = -1.67163;
  const double w_p2x = -0.726050;
  const double wa_m2x = -0.0944815;
  const double wa_p2x = 2.008660;
  const double wz_m2x = -1.108270;
  const double wz_p2x = -0.739106;
  //const double wzz_m2x = -0.743778;
  //const double wzz_p2x = -0.363900;
  const double As_m2x = 0.600062;
  const double As_p2x = 1.047250;
  const double h_m2x = 0.592508;
  const double h_p2x = 0.745060;
  const double Omm_m2x = 0.213166;
  const double Omm_p2x = 0.400481;
  const double sig_m2x = 0.562842;
  const double sig_p2x = 0.865596;
  const double sigxx_m2x = 0.;
  const double sigxx_p2x = 0.;
  */
  const double omb_m2x = 0.0201515;
  const double omb_p2x = 0.0256071;
  const double omc_m2x = 0.0923604;
  const double omc_p2x = 0.121335;
  const double n_m2x = 0.924424;
  const double n_p2x = 1.022070;
  const double w_m2x = -1.31900;
  const double w_p2x = -0.884513;
  const double wa_m2x = 0.2090045;
  const double wa_p2x = 1.451430;
  const double wz_m2x = -0.984745;
  const double wz_p2x = -0.791281;
  //const double wzz_m2x = -0.743778;
  //const double wzz_p2x = -0.363900;
  const double As_m2x = 0.682938;
  const double As_p2x = 0.978227;
  const double h_m2x = 0.628916;
  const double h_p2x = 0.704714;
  const double Omm_m2x = 0.230481;
  const double Omm_p2x = 0.326418;
  const double sig_m2x = 0.692005;
  const double sig_p2x = 0.806707;
  const double sigxx_m2x = 0.;
  const double sigxx_p2x = 0.;
#else
  const double omb_m2x = 0.0196275;
  const double omb_p2x = 0.0263619;
  const double omc_m2x = 0.0917873;
  const double omc_p2x = 0.118932;
  const double n_m2x = 0.912716;
  const double n_p2x = 1.05649;
  const double w_m2x = -1.05071;
  const double w_p2x = -0.755120;
  const double As_m2x = 0.597036;
  const double As_p2x = 0.988167;
  const double h_m2x = 0.653634;
  const double h_p2x = 0.762706;
  const double Omm_m2x = 0.212878;
  const double Omm_p2x = 0.303393;
  const double sig_m2x = 0.689439;
  const double sig_p2x = 0.900127;
  const double sigxx_m2x = 0.;
  const double sigxx_p2x = 0.;
#endif
#endif

  dbg<<"Fully marginalized best fits:\n";
  dbg<<"omega_b = "<<omb_fit<<"  68%: ("<<omb_m1<<","<<omb_p1<<")  95%: ("<<omb_m2<<","<<omb_p2<<")\n";
  dbg<<"     absolute: +"<<omb_p2-omb_fit<<" - "<<omb_fit-omb_m2<<endl;
  dbg<<"     relative: +"<<omb_p2-omb_p2x<<" -"<<omb_m2x-omb_m2<<endl;
  dbg<<"omega_c = "<<omc_fit<<"  68%: ("<<omc_m1<<","<<omc_p1<<")  95%: ("<<omc_m2<<","<<omc_p2<<")\n";
  dbg<<"     absolute +"<<omc_p2-omc_fit<<" -"<<omc_fit-omc_m2<<endl;
  dbg<<"     relative +"<<omc_p2-omc_p2x<<" -"<<omc_m2x-omc_m2<<endl;
  dbg<<"n = "<<n_fit<<"  68%: ("<<n_m1<<","<<n_p1<<")  95%: ("<<n_m2<<","<<n_p2<<")\n";
  dbg<<"     absolute +"<<n_p2-n_fit<<" -"<<n_fit-n_m2<<endl;
  dbg<<"     relative +"<<n_p2-n_p2x<<" -"<<n_m2x-n_m2<<endl;
#ifndef WIsM1
  dbg<<"w = "<<w_fit<<"  68%: ("<<w_m1<<","<<w_p1<<")  95%: ("<<w_m2<<","<<w_p2<<")\n";
  dbg<<"     absolute +"<<w_p2-w_fit<<" -"<<w_fit-w_m2<<endl;
  dbg<<"     relative +"<<w_p2-w_p2x<<" -"<<w_m2x-w_m2<<endl;
#endif
#ifdef WWA
  dbg<<"wa = "<<wa_fit<<"  68%: ("<<wa_m1<<","<<wa_p1<<")  95%: ("<<wa_m2<<","<<wa_p2<<")\n";
  dbg<<"     absolute +"<<wa_p2-wa_fit<<" -"<<wa_fit-wa_m2<<endl;
  dbg<<"     relative +"<<wa_p2-wa_p2x<<" -"<<wa_m2x-wa_m2<<endl;
  dbg<<"w(z="<<ZPIVOT<<") = "<<wz_fit<<"  68%: ("<<wz_m1<<","<<wz_p1<<")  95%: ("<<wz_m2<<","<<wz_p2<<")\n";
  dbg<<"     absolute +"<<wz_p2-wz_fit<<" -"<<wz_fit-wz_m2<<endl;
  dbg<<"     relative +"<<wz_p2-wz_p2x<<" -"<<wz_m2x-wz_m2<<endl;
  //dbg<<"w(z="<<wzz_z<<") = "<<wzz_fit<<"  68%: ("<<wzz_m1<<","<<wzz_p1<<")  95%: ("<<wzz_m2<<","<<wzz_p2<<")\n";
  //dbg<<"     absolute +"<<wzz_p2-wzz_fit<<" -"<<wzz_fit-wzz_m2<<endl;
  //dbg<<"     relative +"<<wzz_p2-wzz_p2x<<" -"<<wzz_m2x-wzz_m2<<endl;
#endif
  dbg<<"A_s = "<<As_fit<<"  68%: ("<<As_m1<<","<<As_p1<<")  95%: ("<<As_m2<<","<<As_p2<<")\n";
  dbg<<"     absolute +"<<As_p2-As_fit<<" -"<<As_fit-As_m2<<endl;
  dbg<<"     relative +"<<As_p2-As_p2x<<" -"<<As_m2x-As_m2<<endl;
  dbg<<"h = "<<h_fit<<"  68%: ("<<h_m1<<","<<h_p1<<")  95%: ("<<h_m2<<","<<h_p2<<")\n";
  dbg<<"     absolute +"<<h_p2-h_fit<<" -"<<h_fit-h_m2<<endl;
  dbg<<"     relative +"<<h_p2-h_p2x<<" -"<<h_m2x-h_m2<<endl;
  dbg<<"Omega_m = "<<Omm_fit<<"  68%: ("<<Omm_m1<<","<<Omm_p1<<")  95%: ("<<Omm_m2<<","<<Omm_p2<<")\n";
  dbg<<"     absolute +"<<Omm_p2-Omm_fit<<" -"<<Omm_fit-Omm_m2<<endl;
  dbg<<"     relative +"<<Omm_p2-Omm_p2x<<" -"<<Omm_m2x-Omm_m2<<endl;
  dbg<<"sigma_8 = "<<sig_fit<<"  68%: ("<<sig_m1<<","<<sig_p1<<")  95%: ("<<sig_m2<<","<<sig_p2<<")\n";
  dbg<<"     absolute +"<<sig_p2-sig_fit<<" -"<<sig_fit-sig_m2<<endl;
  dbg<<"     relative +"<<sig_p2-sig_p2x<<" -"<<sig_m2x-sig_m2<<endl;
  dbg<<"sigma_8(Omega_m/0.3)^1/2 = "<<sigxx_fit<<"  68%: ("<<sigxx_m1<<","<<sigxx_p1<<")  95%: ("<<sigxx_m2<<","<<sigxx_p2<<")\n";
  dbg<<"     absolute +"<<sigxx_p2-sigxx_fit<<" -"<<sigxx_fit-sigxx_m2<<endl;
  dbg<<"     relative +"<<sigxx_p2-sigxx_p2x<<" -"<<sigxx_m2x-sigxx_m2<<endl;
#endif

  return 0;
}
